#### Relationship between phthalates exposures and hyperuricemia in U.S. general population ####
#探讨10种邻苯二甲酸酯代谢物与高尿酸血症的关系
# 2023年3月30日15:05:33 里面用到了 Bayesian kernel machine regression 不了解啊不要浪费时间了 去找找合适的23年文章 我能搞懂的
# 2023年5月24日08:52:15 进行筛选数据
# 根据上次文章的经验 还是不要搞的很复杂的excel 大部分完成还是在里面完成
# 方法 本研究为横断面研究。 研究对象来自2007至2016年NHANES数据库中20岁以上的6865名研究对象。 所有研
#究对象均有10种邻苯二甲酸酯代谢产物(MECPP、 MnBP、 MEHHP、 MEOHP、 MiBP、 cx-MiNP、 MCOP、
# MCPP、 MEP、 MBzP)、 高尿酸血症和协变量的完整数据。 采用多变量logistics回归、 限制性立方样条(RCS)模型
# 和贝叶斯核机器回归(BKMR)模型评估邻苯二甲酸酯代谢物与高尿酸血症之间的单一、 非线性和混合关系。
# 作为补充， 我们还评估了邻苯二甲酸酯代谢物与血清尿酸(SUA)水平之间的关系

# 2007—2016年NHANES共有50,588名参与者， 本研究纳入了26,355名
#20岁以上且有尿酸数据的符合条件的参与者。 由于仅在其中1 / 3的参
#与者中测定了邻苯二甲酸酯代谢物， 因此有8,493名参与者的暴露数据。
#我们排除了1,628名没有协变量数据的参与者。 有可用数据的6,865名参
#与者被纳入统计学分析。 图1显示了参与者选择的流程图
# 2023年7月5日10:07:53 重新梳理下 这篇文章主要是研究了如上面所说的10种邻苯二甲酸酯代谢物(10pm (简写后续识别))与高尿酸血症(hyperuricemia)的关系 
# 还研究了serum uric acid (SUA) 血清尿酸与10pm关系 
# 
library(tidyverse)
library(gtsummary)
library(haven)
library(skimr)
library(corrplot)
library(bkmr)
library(forestplot)
library(rms)

setwd("G:/BaiduNetdiskDownload")
# 2023年6月3日14:14:53根据补充材料 因为自己筛选的变量和补充的部分不一样 重新核对下
#读取2007-2016年数据 
# demo 数据 基础年龄 性别 教育程度等

demo_e <- read_xpt("NHANES/2007-2008/Demographics/demo_e.xpt")
demo_f <- read_xpt("NHANES/2009-2010/Demographics/demo_f.xpt")
demo_g <- read_xpt("NHANES/2011-2012/Demographics/demo_g.xpt")
demo_h <- read_xpt("NHANES/2013-2014/Demographics/demo_h.xpt")
demo_i <- read_xpt("NHANES/2015-2016/Demographics/demo_i.xpt")

demo_all <- dplyr::bind_rows(list(demo_e,demo_f,demo_g,demo_h,demo_i))

# BMI 
bmi_e <- read_xpt("NHANES/2007-2008/Examination/bmx_e.xpt")
bmi_f <- read_xpt("NHANES/2009-2010/Examination/bmx_f.xpt")
bmi_g <- read_xpt("NHANES/2011-2012/Examination/bmx_g.xpt")
bmi_h <- read_xpt("NHANES/2013-2014/Examination/bmx_h.xpt")
bmi_i <- read_xpt("NHANES/2015-2016/Examination/bmx_i.xpt")



bmi_all <- dplyr::bind_rows(list(bmi_e,bmi_f,bmi_g,bmi_h,bmi_i))


# bpx_all <- dplyr::bind_rows(list(bpx_e,bpx_f,bpx_g,bpx_h,bpx_i))

# 实验室数据 尿酸 BIOPRO	lbdsuasi	Uric acid (umol/L)	尿酸 (umol/L)

# 饮酒
alq.e <- read_xpt("NHANES/2007-2008/Questionnaire/alq_e.xpt")
alq.f <- read_xpt("NHANES/2009-2010/Questionnaire/alq_f.xpt")
alq.g <- read_xpt("NHANES/2011-2012/Questionnaire/alq_g.xpt")
alq.h <- read_xpt("NHANES/2013-2014/Questionnaire/alq_h.xpt")
alq.i <- read_xpt("NHANES/2015-2016/Questionnaire/alq_i.xpt")

alq_all <- dplyr::bind_rows(list(alq.e,alq.f,alq.g,alq.h,alq.i))

# 血压问卷
bpq.e <- read_xpt("NHANES/2007-2008/Questionnaire/bpq_e.xpt")
bpq.f <- read_xpt("NHANES/2009-2010/Questionnaire/bpq_f.xpt")
bpq.g <- read_xpt("NHANES/2011-2012/Questionnaire/bpq_g.xpt")
bpq.h <- read_xpt("NHANES/2013-2014/Questionnaire/bpq_h.xpt")
bpq.i <- read_xpt("NHANES/2015-2016/Questionnaire/bpq_i.xpt")

bpq_all <- dplyr::bind_rows(list(bpq.e,bpq.f,bpq.g,bpq.h,bpq.i))

# 活动
paq_e <- read_xpt("NHANES/2007-2008/Questionnaire/paq_e.xpt")
paq_f <- read_xpt("NHANES/2009-2010/Questionnaire/paq_f.xpt")
paq_g <- read_xpt("NHANES/2011-2012/Questionnaire/paq_g.xpt")
paq_h <- read_xpt("NHANES/2013-2014/Questionnaire/paq_h.xpt")
paq_i <- read_xpt("NHANES/2015-2016/Questionnaire/paq_i.xpt")
paq_all <- dplyr::bind_rows(list(paq_e,paq_f,paq_g,paq_h,paq_i))

# CKD 慢性肾脏病
kiq_u_e <- read_xpt("NHANES/2007-2008/Questionnaire/kiq_u_e.xpt")
kiq_u_f <- read_xpt("NHANES/2009-2010/Questionnaire/kiq_u_f.xpt")
kiq_u_g <- read_xpt("NHANES/2011-2012/Questionnaire/kiq_u_g.xpt")
kiq_u_h <- read_xpt("NHANES/2013-2014/Questionnaire/kiq_u_h.xpt")
kiq_u_i <- read_xpt("NHANES/2015-2016/Questionnaire/kiq_u_i.xpt")
kiq_u_all <- dplyr::bind_rows(list(kiq_u_e,kiq_u_f,kiq_u_g,kiq_u_h,kiq_u_i))


#糖尿病 Diabetes
diq_e <- read_xpt("NHANES/2007-2008/Questionnaire/diq_e.xpt")
diq_f <- read_xpt("NHANES/2009-2010/Questionnaire/diq_f.xpt")
diq_g <- read_xpt("NHANES/2011-2012/Questionnaire/diq_g.xpt")
diq_h <- read_xpt("NHANES/2013-2014/Questionnaire/diq_h.xpt")
diq_i <- read_xpt("NHANES/2015-2016/Questionnaire/diq_i.xpt")
diq_all <- dplyr::bind_rows(list(diq_e,diq_f,diq_g,diq_h,diq_i))

biopro_e <- read_xpt("NHANES/2007-2008/Laboratory/biopro_e.xpt")
biopro_f <- read_xpt("NHANES/2009-2010/Laboratory/biopro_f.xpt")
biopro_g <- read_xpt("NHANES/2011-2012/Laboratory/biopro_g.xpt")
biopro_h <- read_xpt("NHANES/2013-2014/Laboratory/biopro_h.xpt")
biopro_i <- read_xpt("NHANES/2015-2016/Laboratory/biopro_i.xpt")
biopro_all <- dplyr::bind_rows(list(biopro_e,biopro_f,biopro_g,biopro_h,biopro_i))

# 尿肌酐
alb_cr_e <- read_xpt("NHANES/2007-2008/Laboratory/alb_cr_e.xpt")
alb_cr_f <- read_xpt("NHANES/2009-2010/Laboratory/alb_cr_f.xpt")
alb_cr_g <- read_xpt("NHANES/2011-2012/Laboratory/alb_cr_g.xpt")
alb_cr_h <- read_xpt("NHANES/2013-2014/Laboratory/alb_cr_h.xpt")
alb_cr_i <- read_xpt("NHANES/2015-2016/Laboratory/alb_cr_i.xpt")
alb_cr_all <- dplyr::bind_rows(list(alb_cr_e,alb_cr_f,alb_cr_g,alb_cr_h,alb_cr_i))

# 血清可替宁 多数据模块拿取
cotnal_e <- read_xpt("NHANES/2007-2008/Laboratory/cotnal_e.xpt")
cotnal_f <- read_xpt("NHANES/2009-2010/Laboratory/cotnal_f.xpt")
cotnal_g <- read_xpt("NHANES/2011-2012/Laboratory/cotnal_g.xpt")
cot_h <- read_xpt("NHANES/2013-2014/Laboratory/cot_h.xpt")
cot_i <- read_xpt("NHANES/2015-2016/Laboratory/cot_i.xpt")
cot_all <- dplyr::bind_rows(list(cotnal_e,cotnal_f,cotnal_g,cot_h,cot_i)) 



# 实验室数据 https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/PHTHTE_J.htm 
# 关于邻苯二甲酸盐 

phthte_e <- read_xpt("NHANES/2007-2008/Laboratory/phthte_e.xpt")
phthte_f <- read_xpt("NHANES/2009-2010/Laboratory/phthte_f.xpt")
phthte_g <- read_xpt("NHANES/2011-2012/Laboratory/phthte_g.xpt")
phthte_h <- read_xpt("NHANES/2013-2014/Laboratory/phthte_h.xpt")
phthte_i <- read_xpt("NHANES/2015-2016/Laboratory/phthte_i.xpt")

phthte_all <- dplyr::bind_rows(list(phthte_e,phthte_f,phthte_g,phthte_h,phthte_i))

output <- plyr::join_all(list(demo_all, biopro_all,phthte_all,cot_all,alq_all,paq_all,kiq_u_all,diq_all,bmi_all,bpq_all,alb_cr_all),
                         by='SEQN',type='full')

# 50588    
dim(output)
# 排除20岁以下的
outputagelower <- output[which(output$RIDAGEYR<20),]
dim(outputagelower) # 21387    
# 不小于20岁的人群  29201    
outputagelarge <- output[which(output$RIDAGEYR>=20),]
# 2846 serum uric acid  血清尿酸  (SUA)
outnaacid <- output[which(is.na(outputagelarge$LBDSUASI)),]
#26355
outputagelargenotnaacid <- outputagelarge[which(!is.na(outputagelarge$LBDSUASI)),]
dim(outputagelargenotnaacid)

# 删除不存在邻苯二甲酸盐 衍生物数据的数据
# All participants had complete data on ten phthalate metabolites (MECPP, MnBP, MEHHP, MEOHP, MiBP, cx-MiNP, MCOP,
# MCPP, MEP, MBzP), hyperuricemia, and covariates

#MECPP ->   mono-2-ethyl-5-carboxypentyl phthalate -> URXECP
#MnBP -> mono-n-butyl phthalate ->  URXMBP
#MEHHP ->  mono(2-ethyl-5-hydroxyhexyl) phthalate (MEHHP) ->  URXMHH
#MEOHP ->  mono(2-ethyl-5-oxohexyl) phthalate (MEOHP) ->  URXMOH
#MiBP -> mono-isobutyl phthalate (MiBP) ->   URXMIB
# cx-MiNP -> mono(carboxynonyl) phthalate ->   URXCNP
#MCOP -> monocarboxyoctyl phthalate ->  URXCOP 
#MCPP -> mono(3-carboxypropyl) phthalate (MCPP) ->  URXMC1
#MEP -> mono ethyl phthalate (MEP) -> URXMEP
#MBzP ->mono benzyl phthalate -> URXMZP

outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXECP)),]
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXMBP)),]
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXMHH)),]
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXMOH)),]
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXMIB)),]
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXCNP)),]
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXCOP)),]
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXMC1)),]
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXMEP)),]
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXMZP)),]
dim(outputagelargenotnaacid) # [1] 8493  132 

# 2023年5月24日16:12:59 对得上 明天搞剩下的 

# blood pressure measures
# body measures 
# serum cotinine
# albumin creatinine 
# oral glucose tolerance test
# plasma fasting glucose
# glycohemoglobin measures
# alcohol use
# blood pressure & cholesterol
# diabetes
# physical activity
# 删除了缺少BMI 血清可替宁  等数据
# BMI -> BMXBMI  Examination\\bmx_j.xpt模块
# serum cotinine -> LBXCOT    Laboratory/cot_e 这个放一边 是两种不同的组合 COT -> LBXCOT COTNAL -> LBXCOT
# Family PIR  -> INDFMPIR 
# Drinking status -> ALQ101 
# education level  -> DMDEDUC2略过!
# physical activity -> PAQ650 PAQ650',' ? 哪个为准呢 略过!
# urine creatinine  -> URXCRS   ALB 模块 实验室数据
# hypertension -> BPQ020 -- Questionnaire/bpq_d
# 2023年5月26日10:43:05 这里部分没有仔细说明 部分实在找不出对应的数量 部分是对的起来的 进行后续操作吧
# 测试EDU 教育程度
# outnaaedu  <- outputagelargenotnaacid[which(is.na(outputagelargenotnaacid$DMDEDUC2)),]
# # 测试bmi
# outnaabmi <- outputagelargenotnaacid[which(is.na(outputagelargenotnaacid$BMXBMI)),]
# # 测试pir 贫困
# outnaapir <- outputagelargenotnaacid[which(is.na(outputagelargenotnaacid$INDFMPIR)),]
# 饮酒 就是只能这样搞 alqy 那个变量无法测试 最小数量也要728 Had at least 12 alcohol drinks per year  alq101 原文附表 提示
outnadrink <- which(is.na(outputagelargenotnaacid$ALQ101))
serumcotininena <- which(is.na(outputagelargenotnaacid$LBXCOT))
DMDEDUC2na <- which(is.na(outputagelargenotnaacid$DMDEDUC2))

length(outnadrink) # 728 对不起来
length(serumcotininena) # 2 对的起来
length(DMDEDUC2na) # 0 对不的起来 2023年6月3日16:17:44 今天核对到这里
# table(outalqtemp$DMDEDUC2)
#table(outputagelargenotnaacid$LBXCOT)
#ctinine <- which(is.na(outputagelargenotnaacid$LBXCOT))
#length(ctinine)

# colnames(outputagelargenotnaacid)

# outnaabpq <- outputagelargenotnaacid[which(is.na(outputagelargenotnaacid$BPQ020)),]

# dim(outnaabmi) # 84 BMI 准确 
# dim(serumcotinine) # 2  serumcotinine 准确
# dim(outnaapir) # 757 pir 准确
# dim(outnaaedu) # 0 不准确 差7 
# dim(outnadrink) # 603 不准确 差118
# dim(outnaabpq) # 0
# table(outputagelargenotnaacid$DMDEDUC2)

# 直接删除了 部分数据复现不了 能删除多少算多少  
# 2023年5月27日13:52:22 目前只是删除了BMI  pir edu(0) drink 
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$BMXBMI)),]
dim(outputagelargenotnaacid)
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$INDFMPIR)),]
dim(outputagelargenotnaacid)
outputagelargenotnaacid <- outputagelargenotnaacid[which(!is.na(outputagelargenotnaacid$URXECP)),]
dim(outputagelargenotnaacid)
outputagelargenotnaacid <- outputagelargenotnaacid[-outnadrink,]
dim(outputagelargenotnaacid)
# 最终数据 outputagelargenotnaacid 6991 
outputdata <- outputagelargenotnaacid
# tabl1  年龄 性别 BMI  种族 教育程度 pir 娱乐活动 吸烟状态 饮酒状态 高血压(hypertension) 糖尿病 (diabetes) CKD

# 娱乐活动 PAQ605 - vigorous 剧烈的工作活动 
# 娱乐活动 PAQ620 - moderate 缓和
# 娱乐活动 剩下就是 - no 不运动 官网没找到变量 不运动的

# 饮酒状态 alq110  1 YES 2 NO

# 吸烟状态  LBDCOTLC 0 吸烟 1 不吸烟  serum cotinine, above 10 ng/ml 大于是吸烟 小于等于 不吸烟
# 高血压 Hypertension  bpq_all ->  BPQ020 1 是 2 不是

# 糖尿病  Diabetes 糖尿病前期（pre-DM）到正常血糖调节（ normal glucose regulation ->NGR）和糖尿病（DM)

# diq160 (pre-DM)  {Have you/Has SP} ever been told by a doctor or other health professional that {you have/SP has} any of the following:prediabetes, impaired fasting glucose, impaired glucose tolerance, borderline diabetes or that {your/her/his} blood sugar is higher than normal but not high enough to be called diabetes or sugar diabetes?
# Diabetes      DIQ010 1 yes
# 剩下 ngr 
#CKD KIQ_U  KIQ022 1 YES 2 NO   Questionnaire 
# urine creatinine  -> URXCRS
# Hyperuricemia LBDSUASI
#was defined as dichotomous with SUA levels ≥416 μmol/L (7.0 mg/dL)
#and ≥ 357 μmol/L (6.0 mg/dL) in male and female 
# 设置数据
# 性别
outputdata$Sex <- ifelse(outputdata$RIAGENDR==1,'male','female')
#bmi
outputdata$BMI <- ifelse(outputdata$BMXBMI<=24.9,'≤ 24.9', '> 24.9')
# 种族
race <- recode_factor(outputdata$RIDRETH1, 
                      `1` = 'Mexican American',
                      `2` = 'Other Hispanic',
                      `3` = 'White',
                      `4` = 'Black',
                      `5` = 'Other Race'
                      
)
outputdata$race <- race
# 教育程度
edu <- recode_factor(outputdata$DMDEDUC2, 
                     `1` = 'Less than high school',
                     `2` = 'Less than high school',
                     `3` = 'High school graduate/GED or equivalent',
                     `4` = 'College or above',
                     `5` = 'College or above'
                     
)
outputdata$edu <- edu
# 贫困程度
outputdata$FPL <- ifelse(outputdata$INDFMPIR<1.0,'< 1.0',
                         ifelse((outputdata$INDFMPIR>=1.0&outputdata$INDFMPIR<2.0) , "1.0–2.0" ,
                                ifelse((outputdata$INDFMPIR>=2.0&outputdata$INDFMPIR<4.0) , "2.0–4.0" ,'> 4.0')))

# 活动
# 娱乐活动 PAQ605 - vigorous 剧烈的工作活动 
# 娱乐活动 PAQ620 - moderate 缓和
# 娱乐活动 剩下就是 - no 不运动 官网没找到变量 不运动的

outputdata$activity <- ifelse(outputdata$PAQ605==1,'vigorous',
                         ifelse((outputdata$PAQ620==1) , "moderate" ,'no'))

# 吸烟 1 是 
outputdata$smoke <- ifelse(outputdata$LBXCOT>10,'yes','no')
# 饮酒
outputdata$drink <- ifelse(outputdata$ALQ101==1,'yes','no')
#Hypertension
outputdata$Hypertension <- ifelse(outputdata$BPQ020==1,'yes','no')

# Diabetes


# 把na 的补充为 no
outputdata[which(is.na(outputdata$smoke)),'smoke'] <-'no'

outputdata$Diabetes <- ifelse(outputdata$DIQ160==1,'Pre-diabetes',
                              ifelse(outputdata$DIQ010==1 , "Diabetes" ,'NGR'))

table(outputdata$Diabetes) # 变量找的不对  后期修复
table(outputdata$DIQ010)

# ckd 
outputdata$CKD <- ifelse(outputdata$KIQ022==1,'yes','no')
outputdata$urinecreatinine <- outputdata$URXCRS

outputdata$Hyperuricemia <- ifelse((outputdata$Sex=='male' & outputdata$LBDSUASI>=416.0),'Yes',
                                   ifelse((outputdata$Sex=='male' & outputdata$LBDSUASI<416.0),'No',
                                   ifelse((outputdata$Sex=='female' & outputdata$LBDSUASI>=357.0),'Yes',
                                   ifelse((outputdata$Sex=='female' & outputdata$LBDSUASI<357.0),'No',""))))




outputdata$Hyperuricemiaint <- ifelse(outputdata$Hyperuricemia=='Yes',1,0)

outputdatasin <-outputdata[,c('SEQN', 'RIDAGEYR', 'Sex', 'BMI','race','edu','FPL', 'activity','smoke','drink','Hypertension','Diabetes','CKD',
                              'urinecreatinine','URXECP','URXMBP','URXMHH','URXMOH','URXMIB','URXCNP','URXCOP','URXMC1','URXMEP','URXMZP','LBDSUASI','Hyperuricemia','Hyperuricemiaint')]
# colnames(outputdata)
# Spearman斯皮尔曼相关性 使用数据
#unique(outputdatasin$SEQN)
# outputdatasin[outputdatasin$SEQN==47083,]


outputdatasinSpearman <-outputdata[,c('URXECP', 'URXMBP', 'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP')]
outputdatasinSpearmanLBDSUASI <-outputdata[,c('URXECP', 'URXMBP', 'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP','LBDSUASI')]
outputdatasinSpearmanLBDSUASIsin <- outputdatasinSpearmanLBDSUASI[c(1:600),]
colnames(outputdatasinSpearman) <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','cx-MiNP','MCOP','MCPP','MEP','MBzP')

outputdatasin$Sex<-as.factor(outputdatasin$Sex)
outputdatasin$BMI<-as.factor(outputdatasin$BMI)
outputdatasin$race<-as.factor(outputdatasin$race)
outputdatasin$edu<-as.factor(outputdatasin$edu)
outputdatasin$FPL<-as.factor(outputdatasin$FPL)
outputdatasin$activity<-as.factor(outputdatasin$activity)
outputdatasin$smoke<-as.factor(outputdatasin$smoke)
outputdatasin$drink<-as.factor(outputdatasin$drink)
outputdatasin$Hypertension<-as.factor(outputdatasin$Hypertension)
outputdatasin$Diabetes<-as.factor(outputdatasin$Diabetes)
outputdatasin$CKD<-as.factor(outputdatasin$CKD)

# skim(outputdatasin)

# tbl_summary(data = outputdatasin, missing = 'no' ,by =Hyperuricemia )%>%add_overall()

###############################
# https://www.zhihu.com/tardis/zm/art/109826800?source_id=1005
##Pearson、Spearman、Kendall 相关
mtcars1 <-data(mtcars)

#标准化不影响相关系数计算值，但可以让数据服从均值 0，标准差 1 的等方差结构
mtcars <- scale(mtcars)

# #协方差计算，cov()
# cov_pearson <- cov(mtcars, method = 'pearson')
# cov_pearson
# 
# cov_spearman <- cov(mtcars, method = 'spearman')
# cov_spearman
# 
# cov_kendall <- cov(mtcars, method = 'kendall')
# cov_kendall

#相关系数计算，cor()
# cor_pearson <- cor(mtcars, method = 'pearson')
# cor_pearson

cor_spearman <- cor(mtcars, method = 'spearman')
colnames(mtcars)
cor_spearmanMy <- cor(outputdatasinSpearman, method = 'spearman')
colnames(outputdatasin)
cor_spearman



#相关图，例如

# #应该就是这个画的 还需要研究参数问题
# corrplot(cor_pearson, method = 'pie', type = 'lower',number.cex = 0.8, diag = FALSE, tl.cex = 1.2,add = F)
# # circle", "square", "ellipse", "number", "shade", "color", "pie"
# corrplot(cor_spearman, add = FALSE, type = 'upper', method = 'circle', diag = FALSE, tl.pos = 'n', cl.pos = 'n')
# corrplot(cor_spearman, add = FALSE, type = 'lower', method = 'square', diag = FALSE, tl.pos = 'n', cl.pos = 'n')
# corrplot(cor_spearman, add = FALSE, type = 'upper', method = 'ellipse', diag = FALSE, tl.pos = 'n', cl.pos = 'n')
# corrplot(cor_spearman, add = FALSE, type = 'upper', method = 'number', diag = FALSE, tl.pos = 'n', cl.pos = 'n')
# corrplot(cor_spearman, add = FALSE, type = 'upper', method = 'shade', diag = FALSE, tl.pos = 'n', cl.pos = 'n')
# corrplot(cor_spearman, add = FALSE, type = 'upper', method = 'color', diag = FALSE, tl.pos = 'n', cl.pos = 'n')
# corrplot.mixed(cor_spearman,  lower = 'pie', upper = 'number')
# corrplot.mixed(cor_spearman,  lower = 'pie', upper = 'number',order = 'hclust',add = F)
# corrplot.mixed(cor_spearman,  lower = 'pie', upper = 'number',order = 'hclust')
# corrplot.mixed(cor_spearman,  lower = 'pie', upper = 'number',order = 'hclust')
# corrplot.mixed(cor_spearman,  lower = 'pie', upper = 'number',order = 'hclust')
# corrplot.mixed(cor_spearman,  lower = 'pie', upper = 'number',order = 'hclust')
# library(ggplot2)
# # 图结束
corrplot.mixed(cor_spearmanMy, order = 'hclust',tl.pos = 'lt',upper = 'number',lower = 'circle')

# p1<-{
#   corrplot.mixed(cor_spearmanMy, order = 'hclust',tl.pos = 'lt',upper = 'number',lower = 'circle')
#   recordPlot()
# }
# 
# ggsave(filename = 'aaa.png',plot = replayPlot(p1))

#输出，例如
# write.table(cor_pearson, 'cor_pearson.txt', sep = '\t', col.names = NA, quote = FALSE)

# 2023年5月31日15:20:14 先生成这是个变量的四分位数

 outputdatasinSpearman <-outputdata[,c('SEQN','RIDAGEYR', 'Sex','BMI', 'race', 'edu', 'FPL', 'activity', 'smoke', 'drink', 'Hypertension', 'Diabetes', 'CKD','urinecreatinine', 'URXECP','URXMBP','URXMHH','URXMOH','URXMIB','URXCNP','URXCOP','URXMC1','URXMEP','URXMZP','LBDSUASI','Hyperuricemia','Hyperuricemiaint')]
 colnames(outputdatasinSpearman) <-  c('SEQN', 'Age', 'Sex','BMI', 'race', 'edu', 'FPL', 'activity', 'smoke', 'drink', 'Hypertension', 'Diabetes', 'CKD','urinecreatinine', 'MECPP',  'MnBP', 'MEHHP',  'MEOHP', 'MiBP', 'cxMiNP','MCOP', 'MCPP',   'MEP',   'MBzP','LBDSUASI','Hyperuricemia','Hyperuricemiaint')
 
 # fit1 <-glm(LBDSUASI ~MECPP , data = outputdatasinSpearman)
 # fit2 <-glm(outputdatasinSpearman$LBDSUASI ~. , data = outputdatasinSpearman)
 # tbl_regression(fit1,exponentiate = T)
 # tbl_regression(fit2,exponentiate = T)
# 加载数据至环境中 rms 包环境要求 
ddist <- datadist(outputdatasin) 
options(datadist="ddist")
# outputdatasin

colnames(outputdatasin)
# rms 包下面的glm 函数包含了很多回归策略 也保留了glm 是一个替代函数 下面的位置 4是经过很多测验的
fitMECPP<- rms::Glm(Hyperuricemiaint ~ RIDAGEYR+Sex+BMI+race+edu+FPL+activity+smoke+drink+ Hypertension+Diabetes+CKD+urinecreatinine+rcs(URXECP,4), data = outputdatasin  )
# 测试贾谦rcs 测试代码########################
ddist<-datadist(outputdatasinSpearman)

options(datadist="ddist")



fit <-lrm(Hyperuricemiaint ~rcs(cxMiNP,nk=5)+Age + Sex + race + edu + FPL + activity + 
            smoke + drink + Hypertension + CKD + urinecreatinine + BMI + 
            MECPP + MnBP + MEHHP + MEOHP + MiBP + MCOP + MCPP + 
            MEP + MBzP,data=outputdatasinSpearman,x=FALSE)
pred_OR<-Predict(fit,cxMiNP,ref.zero=TRUE,fun=exp)
par(mar = c(5, 4, 4, 4) + 0.3)
par(xpd=NA)
ylim.bot<-min(pred_OR[,"lower"])
ylim.top<-max(pred_OR[,"upper"])
plot(pred_OR[,1],pred_OR[,"yhat"], 
     xlab = "MiNP",ylab = paste0("OR"),
     type = "l",ylim = c(ylim.bot,ylim.top),
     col="red",lwd=2) 
lines(pred_OR[,1],pred_OR[,"lower"],lty=2,lwd=1.5)
lines(pred_OR[,1],pred_OR[,"upper"],lty=2,lwd=1.5)
lines(x=range(pred_OR[,1]),y=c(1,1),lty=3,col="grey40",lwd=1.3) #
points(1,1,pch=16,cex=1.2)
####################################

# fitMnBP<- rms::Glm(Hyperuricemiaint ~ MCOP+MCPP+MEP+MEP+rcs(MnBP,4), data = outputdatasinSpearman ,x= TRUE )
# fitMnBP2<- rms::Glm(Hyperuricemiaint ~ MCOP+MCPP+MEP+MEP+rcs(MnBP,4), data = outputdatasinSpearman ,x= TRUE )
# fitMECPP<- rms::Glm(LBDSUASI ~ rcs(MECPP,4), data = outputdatasinSpearman ,x= TRUE )
# fitMECPP<- rms::Glm(LBDSUASI ~ rcs(MECPP,4), data = outputdatasinSpearman ,x= TRUE )
# fitMECPP<- rms::Glm(LBDSUASI ~ rcs(MECPP,4), data = outputdatasinSpearman ,x= TRUE )
# fitMECPP<- rms::Glm(LBDSUASI ~ rcs(MECPP,4), data = outputdatasinSpearman ,x= TRUE )
# fitMECPP<- rms::Glm(LBDSUASI ~ rcs(MECPP,4), data = outputdatasinSpearman ,x= TRUE )
# fitMECPP<- rms::Glm(LBDSUASI ~ rcs(MECPP,4), data = outputdatasinSpearman ,x= TRUE )
# fitMECPP<- rms::Glm(LBDSUASI ~ rcs(MECPP,4), data = outputdatasinSpearman ,x= TRUE )
# fitMECPP<- rms::Glm(LBDSUASI ~ rcs(MECPP,4), data = outputdatasinSpearman ,x= TRUE )
fitMECPP
extractAIC(fitMECPP)
anova(fitMECPP)

or1 <- rms::Predict(fitMECPP,drink,type = "predictions",fun = exp,loglog = T )
or2 <- rms::Predict(fitMECPP,drink,type = "predictions" )
or2 <- rms::Predict(fitMECPP,URXECP,type = "predictions" ,fun = exp,loglog = T)
or3 <- rms::Predict(fitMECPP,URXECP,type = "predictions" )
ggplot(or1)
ggplot(or2)
ggplot(or2)
ggplot(or3)


outputdata$MECPP_quantile <- cut(outputdata$URXECP,breaks = quantile(outputdata$URXECP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))
outputdata$MnBP_quantile <- cut(outputdata$URXMBP,breaks = quantile(outputdata$URXMBP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
outputdata$MEHHP_quantile <- cut(outputdata$URXMHH,breaks = quantile(outputdata$URXMHH),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
outputdata$MEOHP_quantile <- cut(outputdata$URXMOH,breaks = quantile(outputdata$URXMOH),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
outputdata$MiBP_quantile <- cut(outputdata$URXMIB,breaks = quantile(outputdata$URXMIB),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
outputdata$cxMiNP_quantile <- cut(outputdata$URXCNP,breaks = quantile(outputdata$URXCNP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
outputdata$MCOP_quantile <- cut(outputdata$URXCOP,breaks = quantile(outputdata$URXCOP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
outputdata$MCPP_quantile <- cut(outputdata$URXMC1,breaks = quantile(outputdata$URXMC1),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
outputdata$MEP_quantile <- cut(outputdata$URXMEP,breaks = quantile(outputdata$URXMEP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
outputdata$MBzP_quantile <- cut(outputdata$URXMZP,breaks = quantile(outputdata$URXMZP),labels = c('Q1', 'Q2', 'Q3', 'Q4')) 





varis<- c('MECPP_quantile','MnBP_quantile','MEHHP_quantile',
          'MEOHP_quantile','MiBP_quantile','cxMiNP_quantile','MCOP_quantile',
          'MCPP_quantile','MEP_quantile','MBzP_quantile')


sinLogistic <- function(x){
  my_list<-list(c(1))
  for (variable in x) {
    len <- length(my_list)
    f<- as.formula(paste0('Hyperuricemiaint ~ ',variable,'+URXCRS'))
    glmMECPP<-glm(f,data = outputdata,family =  binomial )
    my_list[[len+1]] <- glmMECPP
  }
  return(my_list)
}
allSin <-sinLogistic(varis)
# allsinregre <-tbl_stack(tbls = list(allSin[[2]],allSin[[3]],allSin[[4]],allSin[[5]],allSin[[6]],allSin[[7]],allSin[[8]],allSin[[9]],allSin[[10]],allSin[[11]]) )
# 准备森林图数据
my_listF <- list(c(1))

for (variable in c(2:11)) {
  print(variable)
  len <- length(my_listF)
  fit.result<-summary(allSin[[variable]])
  df1<-fit.result$coefficients
  df2<-confint(allSin[[variable]])
  df3<-cbind(df1,df2)
  df4<-data.frame(df3[-1,c(1,4,5,6)])
  df4$Var<-rownames(df4)
  colnames(df4)<-c("OR","Pvalue","OR_1","OR_2","Var")
  df5<-df4[,c(5,1,2,3,4)]
  df5$OR_mean<-df5$OR
  df5$OR<-paste0(round(df5$OR,2),
                 "(",
                 round(df5$OR_1,2),
                 "~",
                 round(df5$OR_2,2),
                 ")")
  df5$Pvalue<-round(df5$Pvalue,3)
  my_listF[[len+1]] <- df5
}

# 多个表进行合并到一个大表显示
my_listF[[3]]
my_listF
allfros <- dplyr::bind_rows(list(my_listF[[2]],my_listF[[3]],my_listF[[4]],my_listF[[5]],my_listF[[6]],my_listF[[7]],my_listF[[8]],my_listF[[9]],my_listF[[10]],my_listF[[11]]))
allfros


write.csv(df5,file = "forestplot_example.csv",
          quote = F,row.names = F)
# 准备结束

fp<-allfros

forestplot(labeltext=as.matrix(fp[,1:3]),
           mean=fp$OR_mean,
           lower=fp$OR_1,
           upper=fp$OR_2,
           zero=0.25,
           boxsize=0.2,
           graph.pos=2)
#  剩下明天搞! 2023-6-2 16:19:42

#1 

#2 




# 表格最后测试复现 好像不是这么回事 看下油管视频啊
# http://www.idata8.com/rpackage/bkmr/kmbayes.html
# https://www.youtube.com/watch?v=8bNFUFOfExU

# set.seed(111)
# fitkm <- kmbayes(y = outputdatasinSpearmanLBDSUASIsin$LBDSUASI, 
#                  Z = outputdatasinSpearmanLBDSUASIsin[,c('URXECP', 'URXMBP', 'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP')], 
#                  iter = 10, verbose = FALSE, 
#                  varsel = TRUE)
# 
# risks.overall <- OverallRiskSummaries(fit = fitkm, y = outputdatasinSpearmanLBDSUASIsin$LBDSUASI, 
#                                       Z = outputdatasinSpearmanLBDSUASIsin[,c('URXECP', 'URXMBP', 'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP')],
#                                       qs = seq(0.25, 0.75, by = 0.05), 
#                                       q.fixed = 0.25, method = "approx")
# risks.overall
# 
# ggplot(risks.overall, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + geom_pointrange()
######################################### 2023年6月2日10:15:12 #####################################
# 根据ai-gpt 给的协变量的bayes 模型 测试
 # 目前 没找到
# Load required packages
# library(bkmr)
# library(randomForest)
# library(glmnet)
# library(HDInterval)
# library(vcd)
# 
# # Load data
# data <- outputdata

#  年龄、 性别、 人种/族群、 贫困状况(PIR)、
#教育水平、 体重指数(BMI))、 生活方式(体力活动、 吸烟状况(血清可替
 #                                 宁)、 饮酒状况)、 病史(糖尿病、 高血压、 慢性肾脏病(CKD))和尿肌酐
#浓度作为协变量
# Set exposure and outcome variables
# 文章中的暴露变量是什么?
# exposure <- data$exposure
# 
# outcome <- data$Hyperuricemiaint
# 
# # Set covariates
# covariates <- data[,c("RIDAGEYR", "Sex", "race","race","FPL","edu","BMI","activity","smoke","drink","Diabetes","Hypertension","CKD","urinecreatinine")]
# 
# covariates <- covariates[c(1:200),]
# outcome <- (outcome[c(1:200)])
# # 暂时搞不定,放一边了 2023年6月2日13:30:17 
# 
# set.seed(111)
# mixture<-as.matrix(outputdatasinSpearman)
# mixture <- mixture[c(1:200),]
# fitkm <- kmbayes(y = outcome, Z = mixture, iter = 5, verbose = FALSE, varsel = TRUE,family = binomial)
# fitkm
# ExtractPIPs(fitkm)
# 
# TracePlot(fit = fitkm, par = "beta")
# TracePlot(fit = fitkm, par = "sigsq.eps")
# TracePlot(fit = fitkm, par = "r", comp = 1)
# 
# risks.overall <- OverallRiskSummaries(fit = fitkm, y = outcome, Z = mixture,
#                                       qs = seq(0.25, 0.75, by = 0.05), 
#                                       q.fixed = 0.5, method = "exact")
# 
# ggplot(risks.overall, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + geom_pointrange()
# 

